function [p, share_t, simShare_t, simSharep_t, simSharep_t_weighted, ...
    dsdp, elas] = equi_p0(p0, w, exp_u_wo_price,...
    ind_endogenous, carrier_t, sim_pcoeff, is_same_firm, setup)

global diff_iter

n_smln = max([size(sim_pcoeff, 2), size(exp_u_wo_price, 2)]);
diff_max = 1;
n_max_iter = 500;
diff_iter = diff_max*ones(n_max_iter, 1);
n_prod_t = size(p0, 1);
% n_carrier = size(carrier_t, 2);
p = p0;
prev_p = p;
count = 0;
% 
% p_ub = 40;
% 
% p_lb = 0.1;


if isscalar(sim_pcoeff) || size(sim_pcoeff, 1)>1
    sim_pcoeff_t = sim_pcoeff;
elseif size(sim_pcoeff, 1)==1
    sim_pcoeff_t = repmat(sim_pcoeff, [n_prod_t 1]);
end

while diff_max > 1e-8 && count<5000
    if setup.use_p0 && count >= 1 % use p from the data
        break;
    end
    
    prev_p = p; % save last round's p       
    
    expu = bsxfun(@times, exp(p.*sim_pcoeff_t), exp_u_wo_price);
    sh_rc = expu./(1+sum(expu, 1));
    
    sh = mean(sh_rc, 2);
    
    dsdp_part1 = diag(mean(sh_rc.*sim_pcoeff_t, 2));
    
    dsdp_part2  = sh_rc*(sh_rc.*sim_pcoeff_t)'.*is_same_firm/n_smln;
    
    p = w + dsdp_part1\dsdp_part2'*(p - w) - dsdp_part1\sh; 
    
    count = count+1;
    diff_max = max(abs(p-prev_p));
    diff_iter(count) = diff_max;
    
%     p = prev_p;
    p = p*0.9 + prev_p*0.1;
end

if count>=5000
    error('price equi not found');
end

p = prev_p;

expu = bsxfun(@times, exp(p.*sim_pcoeff), exp_u_wo_price);
simShare_t = expu./repmat(1+sum(expu, 1), [n_prod_t 1]);

share_t = mean(simShare_t, 2);

simSharep_t = simShare_t.*sim_pcoeff_t;
simSharep_t_weighted = simSharep_t;

simSharep_t_weighted = simSharep_t_weighted ./ n_smln;

dsdp = diag(sum(simSharep_t_weighted,2)) - (simSharep_t_weighted)*simShare_t';


elas = dsdp.*p'./share_t;



% purpose: compute the carrier pricing eqm
% global diff_iter
% 
% n_smln = max([size(sim_pcoeff, 2), size(exp_u_wo_price, 2)]);
% diff_max = 1;
% n_max_iter = 500;
% diff_iter = diff_max*ones(n_max_iter, 1);
% n_prod_t = size(p0, 1);
% n_carrier = size(carrier_t, 2);
% p = p0;
% prev_p = p;
% count = 0;
% 
% p_ub = 40;
% 
% p_lb = 0.1;
% 
% if isscalar(sim_pcoeff) || size(sim_pcoeff, 1)>1
%     sim_pcoeff_t = sim_pcoeff;
% elseif size(sim_pcoeff, 1)==1
%     sim_pcoeff_t = repmat(sim_pcoeff, [n_prod_t 1]);
% end
% 
% while diff_max > 1e-6 && count<5000
%     if setup.use_p0 && count >= 1 % use p from the data
%         break;
%     end
%     
%     prev_p = p; % save last round's p
%     
%     expu = bsxfun(@times, exp(p.*sim_pcoeff), exp_u_wo_price);
%     simShare_t = expu./repmat(1+sum(expu, 1), [n_prod_t 1]);
% 
%     share_t = mean(simShare_t, 2);
%     
%     simSharep_t = simShare_t.*sim_pcoeff_t;
%     simSharep_t_weighted = simSharep_t;
% 
%     simSharep_t_weighted = simSharep_t_weighted ./ n_smln;
% 
%     dsdp = diag(sum(simSharep_t_weighted,2)) - (simSharep_t_weighted)*simShare_t'; 
% 
%     if isfield(setup, 'collusionweightmatrix')
%         if all(ind_endogenous)
%             p = w - (setup.collusionweightmatrix .* dsdp')\share_t;
%         else
%             p(ind_endogenous) = w(ind_endogenous) - ...
%                 (setup.collusionweightmatrix(ind_endogenous,ind_endogenous) .* dsdp(ind_endogenous,ind_endogenous)')\...
%                 (share_t(ind_endogenous)+dsdp(~ind_endogenous, ~ind_endogenous)'*(p(~ind_endogenous)-w(~ind_endogenous)));
%         end
%     else
%         for nc = 1 : n_carrier        
%             Jct_endo = carrier_t(:, nc) & ind_endogenous;
%             Jct_exog = carrier_t(:, nc) & ~ind_endogenous;
%             if any(Jct_endo) && ~any(Jct_exog)
%                 p(Jct_endo) = w(Jct_endo) - dsdp(Jct_endo, Jct_endo)'\share_t(Jct_endo);
%             elseif any(Jct_endo) && any(Jct_exog)
%                 p(Jct_endo) = w(Jct_endo) - dsdp(Jct_endo, Jct_endo)'\(share_t(Jct_endo)+dsdp(Jct_exog, Jct_endo)'*(p(Jct_exog)-w(Jct_exog)));
%             end
%         end
%     end
%     
%     if any(p>p_ub)
%         p(p>p_ub) = rand(sum(p>p_ub), 1).*(p_ub-p_lb)+p_lb;        
%     end
%     
%     if any(p<p_lb)
%         p(p<p_lb) = rand(sum(p<p_lb), 1).*(p_ub-p_lb)+p_lb;
%     end
%     
%     count = count+1;
%     diff_max = max(abs(p-prev_p));
%     diff_iter(count) = diff_max;
%     
%     p = p*0.8 + prev_p*0.2;
% end
% 
% if count>=5000
%     error('price equi not found');
% end
% 
% p = prev_p;
% 
% elas = dsdp.*p'./share_t;

